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Quantum Monte Carlo Calculation of the Binding Energy of Bilayer Graphene 
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We report diffusion quantum Monte Carlo calculations of the interlayer binding energy of bilayer 
graphene. We find the binding energies of the AA- and AB-stacked structures at the equilibrium 
separation to be 11.5(9) and 17.7(9) meV/atom, respectively. The out-of-plane zone-center optical 
phonon frequency predicted by our binding-energy curve is consistent with available experimental 
results. As well as assisting the modeling of interactions between graphene layers, our results will 
facilitate the development of van der Waals exchange-correlation functionals for density functional 
theory calculations. 

PACS numbers: 61.48.Gh, 71.15.Nc, 02.70.Ss 


van der Waals (vdW) interactions play a crucial role 
in a wide range of physical and biological phenomena, 
from the binding of rare-gas solids to the folding of pro¬ 
teins. Significant efforts are therefore being made to de¬ 
velop computational methods that predict vdW contri¬ 
butions to energies of adhesion, particularly for materi¬ 
als such as multilayer graphene. This task has proved to 
be challenging, however, because vdW interactions are 
caused by nonlocal electron correlation effects. Standard 
first-principles approaches such as density functional the¬ 
ory (DFT) with local exchange-correlation functionals 
do not describe vdW interactions accurately. One tech¬ 
nique for including vdW interactions in a first-principles 
framework is to add energies obtained using pairwise in¬ 
teratomic potentials to DFT total energies; this is the 
so-called DFT-D scheme [IHl. The development of vdW 
density functionals (vdW-DFs) that can describe vdW 
interactions in a seamless fashion is another promising 
approach [^-[^. DFT-based random-phase approxima¬ 
tion (RPA) calculations of the correlation energy ilii 
provide a more sophisticated method for treating vdW 
interactions; however, RPA atomization energies are typ¬ 
ically overestimated by up to 15% for solids ll|, Il2| . 
and hence the accuracy of this approach is unclear. 
Symmetry-adapted perturbation theory based on DFT 
allows one to calculate the vdW interactions between 
molecules and hence, by extrapolation, between nanos¬ 
tructures [l3| . Finally, empirical interatomic potentials 
with tails may be used to calculate binding energies 
imi, although such potentials give a qualitatively in¬ 


correct description of the interaction of metallic or tt- 
bonded two-dimensional (2D) materials at large separa¬ 
tion flHl. 


A key test system for methods purporting to describe 
vdW interactions between low-dimensional materials is 
bilayer graphene (BLG). Several theoretical studies have 
used methods based on DFT to calculate the binding 
energy (BE) of BLG. Some of the results are summa¬ 
rized in Table U but there is very little consensus. In 
this work we provide diffusion quantum Monte Garlo 
(DMC) data for the BE of BLG and the atomization 


TABLE I: BE of BLG (both AA- and AB-stacked) obtained 
in recent theoretical studies. The layer separations d quoted 
in the table are the ones used in the calculations, not nec¬ 
essarily the optimized bond length for the given method. 
“SAPT(DFT)” and “DFT-LCAO-00” denote symmetry- 
adapted perturbation theory based on DFT and linear com¬ 
bination of atomic orbitals-orbital occupancy based on DFT, 
respectively. 


Stacking 

Method 

d{A) 

BE (meV/atom) 

AA 

vdW-DF [17] 

3.35 

10.4 

AA 

DFT-D [17] 

3.25 

31.1 

AA 

DMC (pres, wk.) 

3.495 

11.5(9) 

AB 

DFT-LCAO-00 [1^ 

3.1-3.2 

70(5) 

AB 

SAPT(DFT) [1^ 

3.43 

42.5 

AB 

vdW-DF [7] 

3.6 

45.5 

AB 

vdW-DF [17] 

3.35 

29.3 

AB 

DFT-D [17] 

3.25 

50.6 

AB 

DFT-D [2^ 

3.32 

22 

AB 

DMC (pres, wk.) 

3.384 

17.7(9) 


energy of monolayer graphene (MLG), which we have 
extrapolated to the thermodynamic limit. We find the 
DMC BE of BLG to be somewhat less than the BEs pre¬ 
dicted by DFT-D, although the latter vary significantly 
from scheme to scheme. The DMC method is the most 
accurate first-principles technique available for studying 
condensed matter. Our data can therefore be used as a 
benchmark for the development of vdW functionals. 


We have used the variational quantum Monte Carlo 
and DMC methods as implemented in the CASINO code 


2lj to study MLG and BLG. In the former method, 


Monte Carlo integration is used to evaluate expectation 
values with respect to trial many-body wave-function 
forms that may be of arbitrary complexity. In the DMC 
method (23, 


l23l |. a stochastic process governed by the 
Schrodinger equation in imaginary time is simulated to 
project out the ground-state component of the trial wave 
function. Fermionic antisymmetry is maintained by the 
fixed-node approximation, in which the nodal surface is 
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constrained to equal that of the trial wave function [24 1. 
DMC methods have recently been used to study the BE 
of hexagonal boron nitride bilayers 2^ . 

Our many-body trial wave-function form consisted of 
Slater determinants for spin-up and spin-down electrons 
multiplied by a symmetric, positive Jastrow correla¬ 
tion factor exp(J) [i^. The Slater determinants con¬ 
tained Kohn-Sham orbitals that were generated using the 
CASTEP plane-wave DFT code [ 2 ^ within the local den¬ 
sity approximation (LDA). We performed test DMC cal¬ 
culations for 3 X 3 supercells of MLG and AB-stacked 
BLG using Perdew-Burke-Ernzerhof (PBE) orbitals. 
The effect of changing the orbitals on the DMC total en¬ 
ergies (and hence the BE) was statistically insignificant. 

To improve the scaling of our DMC calculations and 
to allow the use of 2D-periodic boundary conditions, 
the orbitals were re-represented in a B-spline (blip) ba¬ 
sis [ 2 ^. The Jastrow exponent J consisted of polyno¬ 
mial and plane-wave expansions in the electron-ion and 
electron-electron distances [ 2 ^. The free parameters in 
the Jastrow factor were optimized by unreweighted vari¬ 
ance minimization [s^, Isij. The DMC energy was ex¬ 
trapolated linearly to zero time step and we verified that 
finite-population errors in our results are negligible ( 3 ^ . 
The fixed-node error is of uncertain magnitude, but it is 
always positive, and should largely cancel when the BE is 
calculated. We used Dirac-Fock pseudopotentials to rep¬ 
resent the C atoms [illsi and fixed the in-plane lattice 
parameter at the experimental value of a = 2.460 A. 

The principal source of uncertainty in our BE results is 
the need to use finite simulation cells subject to periodic 
boundary conditions in DMC calculations for condensed 
matter. Finite-size errors in DMC total energies con¬ 
sist of (i) pseudorandom, oscillatory single-particle finite- 
size errors due to momentum quantization and (ii) sys¬ 
tematic finite-size errors due to the inability to describe 
long-range two-body correlations and the difference be¬ 
tween 1 jr and the 2D Ewald interaction 35, 3^ in a finite 
periodic cell. By dividing the electron-electron interac¬ 
tion energy into a Hartree term (the electrostatic energy 
of the charge density) and an exchange-correlation en¬ 
ergy (the interaction energy of each electron with its ac¬ 
companying exchange-correlation hole) and considering 
the long-range nonoscillatory behavior of the hole pre¬ 
dicted by the RPA, it can be shown that the systematic 
finite-size error in the interaction energy per electron of a 
2D-periodic system is negative and scales asymptotically 
with system size N as [^. The leading-order 

long-range finite-size error in the kinetic energy per elec¬ 
tron behaves in a similar fashion. The finite-size error in 
the atomization energy is therefore positive and scales as 
and the finite-size error in the BE per atom 
must also exhibit the scaling. We also inves¬ 

tigated finite-size errors in the asyniptotic BE using the 
Lifshitz theory of vdW interactions [s^, [s^ with a Dirac 
model of electron dispersion in graphene. To study hnite 


system sizes, we introduced a cutoff wavelength that de¬ 
pended on the cell size and layer separation. However, 
near the equilibrium separation, short-range interactions 
are important and the contribution to the finite-size error 
from the Lifshitz theory is negligible. In order to elim¬ 
inate finite-size effects and obtain the atomization and 
binding energies in the thermodynamic limit, we studied 
simulation cells consisting of arrays of 3 x 3, 4 x 4, and 
6x6 primitive cells for MLG and BLG at the equilib¬ 
rium layer separation and 3x3 and 5x5 cells for BLG 
at nonequilibrium layer separations. We used canonical- 
ensemble twist averaging (40| (i.e., averaging over offsets 
to the grid of k vectors) to reduce the oscillatory single¬ 
particle finite-size errors in the ground-state energies of 
MLG and BLG. To obtain the twist-averaged energy of 
MLG in a simulation cell containing Np primitive cells, 
we performed DMG calculations at twelve random offsets 
kg to the grid of k vectors, then fitted 

E(7Vp,ks) = ij(A^p)-|-6[ALDA(fVp,ks) —Elda(oo)] (1) 

to the DMG energies per atom E{Np,'ks). The model 
function has two fitting parameters: E{Np), which 
is the twist-averaged DMG energy per atom, and b. 
EhDA{Np,'ks) is the DFT-LDA energy per atom of 
MLG obtained using the offset k-point grid correspond¬ 
ing to the supercell used in the DMC calculations, and 
Alda (00) is the DFT-LDA energy per atom obtained us¬ 
ing a fine (50x 50) k-point mesh. Finally, we extrapolated 
our total-energy data to infinite system size by fitting 

E{Np) = E{oo) + (2) 


to the twist-averaged energies per atom, where the ex¬ 
trapolated energy per atom E{oo) and c are fitting pa¬ 
rameters. The atomization energy of MLG is the differ¬ 
ence between the energy of an isolated, spin-polarized C 
atom and the energy per atom of MLG. 

Our DMC atomization energies of MLG as a function 
of system size are plotted in Fig. [T] We find the static- 
nucleus DMC atomization energy to be 7.395(3) eV/atom 
with a Slater-Jastrow trial wave function. This is lower 
than the DMC result of 7.464(10) eV/atom reported in 
Ref. 4l|. Most of this disagreement arises from the use 
of different pseudopotentials in the two works jsij. The 
DFT-PBE phonon zero-point energy (ZPE) of MLG was 
calculated usin g th e method of finite displacements in a 
6x6 supercell [^ and found to be 0.165 eV/atom. The 
ZPE is a correction to be subtracted from the static- 
nucleus atomization energy. In principle, an accurate 
first-principles atomization energy for graphene could be 
used to estimate the BE of graphite by taking the differ¬ 
ence of the experimental atomization energy of graphite 
[7.371(5) eV/atom [i^] and the ZPE-corrected atomiza¬ 
tion energy of MLG. However, the spread of DFT atom¬ 
ization energies resulting from different choices of pseu¬ 
dopotential (of order 40-70 meV/atom [^) implies that 
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FIG. 1: (Color online) Twist-averaged (TA) and non-TA at¬ 
omization energies of MLG against as calculated by 

DMC, where Np is the number of primitive cells in the sim¬ 
ulation supercell. 


FIG. 2: (Color online) Twist-averaged (TA) BLG BE against 
as calculated by DMC, where Np is the number of 
primitive cells in the simulation supercell. The inset shows 
non-twist-averaged BEs. The layer separations are the vdW- 
DF [ 5 ^ equilibrium values of 3.495 and 3.384 A for the AA- 
and AB-stacked structures, respectively. 


first-principles pseudopotential calculations cannot cur¬ 
rently be used to calculate the BE of graphite by this 
approach. 

Despite a great deal of theoretical and experimental 
work, the BE of graphene layers remains poorly under¬ 
stood. The cleavage energy of graphite has been mea¬ 
sured to be 43(5) meV/atom [l^, the BE to be 35(10) 
meV/atom M, and the exfoliation energy to be 52(5) 
meV/atom [^|. More recent experimental work has 
found the cleavage energy to be 31(2) meV/atom (i^ . 
It has been suggested that the latter result may be sub¬ 
stantially underestimated, because the experimental data 
were analyzed using a Lennard-Jones potential, which 
gives qualitatively incorrect interlayer BEs at large sep¬ 
aration 41^. Similar difficulties of interpretation may af¬ 


fect the other experimental results in the literature. The 
results obtained in these works are widely scattered. The 
DMC method has previously been applied to calculate 
the BE of graphite which was found to be 56(5) 
meV/atom, although these calculations were performed 
in relatively small simulation supercells, and finite-size 
effects may limit the accuracy of the results obtained. 

Eor BLG, we restrict our attention to the nonretarded 
regime [d^, in which the BE is simply the difference be¬ 
tween the nonrelativistic total energy per atom in the 
monolayer and the bilayer. We used vdW-DE layer sep¬ 
arations of d = 3.495 A and 3.384 A for the AA- and 
AB-stacked configurations, respectively. In Eig.[2]we plot 
the twist-averaged BEs of AA- and AB-stacked BLG as 
a function of system size. Non-twist-averaged BEs are 
shown in the inset to Eig.[2]and, as expected, show large 
oscillations due to momentum-quantization effects. For 
widely separated graphene layers with nonoverlapping 
charge densities, single-particle finite-size errors cancel 


perfectly when the BE is calculated. However, when the 
layers are closer together, the cancellation is no longer 
perfect. In practice, near the equilibrium separation, 
the single-particle errors in the BE correlate closely with 
the single-particle errors in the total energy of BLG. To 
evaluate the BE in the thermodynamic limit, we twist- 
averaged the BE using Eq. ([T|) with the BE per atom in 
place of E{Np,'ks) and the DFT-LDA total energy per 
atom of BLG in place of kg). We then extrap¬ 

olated the twist-averaged BE to infinite system size using 
Eq. (121). As shown in Fig. [21 the BE of AB-stacked BLG 
is larger than that of AA-stacked BLG, confirming that 
the former is the more stable structure. 

The area of a simulation cell with Np unit cells is 
A = \/iNpa?/2, where a is the lattice parameter of 
graphene. If we define the linear size L of the cell via 
ttL^ = A then we may express the twist-averaged BE 
per atom as Abind(T) = Ebind(oo) -I- where c' is 

—0.31(5) and —0.43(5) eV A®/^ for the AA-stacked and 
AB-stacked geometries, respectively. The BE is reduced 
at small supercell sizes L. The use of a finite supercell 
crudely models the situation where the Gouloinb inter¬ 
action between electrons is screened by a metallic sub¬ 
strate. Hence a metallic substrate is expected to weaken 
the binding of BLG. 

In Fig. [3] we plot the BE of AB-stacked BLG against 
the interlayer separation, as calculated by DFT, DFT-D, 
and DMG. The layer separations we have studied are 
not in the asymptotic regime in which the BE falls off as 
where d is the interlayer separation jsij. We have 
fitted the function Ebind(d) = Aid~^+A^d~^+Ai 2 d~^"^ + 
to our DMG BE data, where the {Ai} are fitting 
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Interlayer distance (A) 

FIG. 3: (Color online) BE curve of AB-stacked BLG as a 
function of interlayer distance calculated using DFT, DFT- 
D, and DMC methods. Our DFT-D calculations used the 
Tkatchenko-Scheffler (TS)_0], Ortmann-Bechstedt-Schmidt 
(OBS) [l^l, and Grimme [5^ vdW corrections. 


parameters, which we find to be A 4 = —2.9 x 10^ meV 
Ag = -2.97x10*5 meVAs, A 12 = 6.18x lO^' meV A^A and 
Aiq = —1.63 X 10® ineV A^®. This function fits the DMC 
data well, with a value of 0.007 per data point. The 
BE found at the minimum of the fitting curve is 17.8(8) 
meV/atom at the equilibrium separation of 3.43(4) A. 
Although the separation that minimizes our fitted BE 
curve for AB-stacked BLG is somewhat larger than the 
separation used in our calculation of the BE reported in 
TablelH the difference between the BEs is not statistically 
significant. The Tkatchenko-Scheffler DFT-D scheme 
shows roughly the same equilibrium separation as DMC, 
but the magnitude of the BE is substantially larger. In 
general, the three DFT-D methods studied [1, M, 
disagree with each other and with DMC. Indeed, the 
magnitude of the BE (if not the shape of the BE curve) 
is best described by the LDA. Our fitted BE curve en¬ 
ables us to calculate the out-of plane zone-center optical 
phonon frequency lozo' of AB-stacked BLG . A com¬ 
parison of ujy.p' frequencies obtained by DFT, DMC, and 
experiments [5^[^ is shown in Table HIl Our DFT-LDA 
frequency is in reasonabl e ag reement with the result (76.8 
cm“5) reported in Ref. The difference between the 
wzo' frequency predicted by our fit to our DMC data and 
the experimental result is negligible [3(7) cm”^] [^ . 


In summary, we have used the DMC method to de¬ 
termine the BE of BLG. Our approach includes a full, 
first-principles treatment of vdW interactions. We have 
found the static-nucleus atomization energy of MLG to 
be 7.395(3) eV/atom, although the uncertainty in this 
result due to the use of nonlocal pseudopotentials may 
be as much as 70 meV/atom [^. We find the BEs of 
AA- and AB-stacked BLG near their equilibrium separa- 


TABLE II: The equilibrium separation do, static-lattice BE at 
equilibrium separation, and out-of-plane zone-center optical- 
phonon frequency uJzo' of AB-stacked BLG obtained by DFT, 
DFT-D, DMG, and experiment. The minimum of the curve 
fitted to the DMC BE data, which is reported in this table, is 
in statistical agreement with the DMC BE obtained using a 
fixed layer separation of 3.384 A, which is reported in Table 


Method 

do (A) 

BE (meV/at.) 

1 uzo’ (cm“^) 

DFT-PBE 

4.40 

1.53 

16 

DFT-LDA 

3.28 

13.38 

84 

DFT-D (TS) 

3.35 

38.03 

111 

DFT-D (OBS) 

3.15 

62.70 

133 

DFT-D (Grimme) 

3.25 

27.08 

95 

DMC (pres, wk.) 

3.43(4) 

17.8(8) 

83(7) 

Exp. 



80(2) 

Exp. 



89 


tions to be 11.5(9) and 17.7(9) meV/atom, respectively. 
Our results indicate that current DFT-D and vdW-DF 
methods significantly overbind 2D materials. 
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